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Abstract 

A numerical technique is described that can efficiently compute solutions in interface prob- 
lems. These are problems with data, such as the coefficients of differential equations, discon- 
tinuous or even singular across one or more interfaces. A prime example of these problems 
are optical waveguides and as such the scheme is applied to Maxwell's equations as they 
are formulated to describe light confinement in Bragg fibers. It is based on standard finite 
differences appropriately modified to take into account all possible discontinuities across the 
waveguide's interfaces due to the change of the refractive index. Second and fourth order 
schemes are described with additional adaptations to handle matrix eigenvalue problems, 
demanding geometries and defects. 
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Light confinement due to cylindrical Bragg reflection instead of total internal reflection was 
first proposed more than three decades ago [l| and gave birth to the so-called Bragg fibers. 
These fibers attract considerable interest because of their ability to guide light in an air 
core; they are essentially dielectric coaxial fibers comprised of alternating circular layers 
with different indices of refraction. The key to making these fibers confine light efficiently, 
i.e., have low absorption loss and a high threshold power for nonlinear effects, is to use 
materials with a high index contrast 0,0, 0, IE!- However, the high index contrast and the 
layered structure that gives these fibers their unique properties also makes them difficult to 
model. 

Mathematically, these problems are called interface problems since their input data (such 
as the coefficients of differential equations, source terms etc.) may be discontinuous or even 
singular across one or several interfaces. The solution to an interface problem, therefore, 
typically is non-smooth or even discontinuous across the interfaces. Interface problems occur 
in many physical applications, particularly for free boundary /moving interface problems, 
such as, the modeling of the Stefan problem of solidification process and crystal growth, 



composite materials, multi-phase flows, cell and bubble deformation, and many others [14 



Several methods have been proposed to study these problems, including asymptotic 
analysis jij], the transfer matrix method Q, finite element methods 0, 0], special func- 
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and Hankel 11 - formalism, and Galerkin numerical methods 12 



tions -Bessel i 

A comparative analysis of the most commonly used methods has also been published [7] 
demonstrating the capabilities and limitations of each method. Among the different nu- 
merical solution methods, the finite difference (FD) method is more attractive due to its 
advantage of simple formulation and numerical implementation and thus will be used here 
to analyze these problem. 



Our approach is based on the immersed interface method 15] (IIM) which has been 
attracting considerable attention due to the many physical problems that can be applied on 



16l Ha [18|, [19|, |2fl, I2JJ22I, |23|, |24 



Galerkin methods 



12 



This approach has two additional advantages over standard 
First of all, the scheme does not need to be modified significantly 
if different boundary conditions are used, thus allowing to calculate all possible solutions 
without any modifications. Methods based upon the Galerkin method typically require a set 
of basis functions that naturally satisfy the boundary conditions, hence the solution must be 
reformulated in a significant way if these change. More importantly, the IIM does not depend 
on any specific functional representation of solutions. Hence, cumbersome integrations or 
finding roots of nontrivial functions, such as Bessel functions (even when asymptotically 
approximated [§[), are avoided. 

The essence of the method is to appropriately modify the correct matrix elements of a 
standard (central) FD scheme so as to take into account all discontinuities across interfaces. 
Starting with a differential equation and under a FD approximation one transforms the 
equation into an algebraic system of the form 



where x is the solution or 



Ax 



Ax = Ax 



for eigenvalue problems, where A is the eigenvalue. The matrix A is comprised of zero 
elements except for the main, upper and lower diagonal (for second order accurate solutions 
or more for more accurate schemes), i.e. 
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This feature gives the FD method an additional advantage since A is tridiagonal. Using 
sparse matrix algebra one can significantly lower computational time, whether the inverse 
of the matrix or the eigenvalues is shout for. The goal is to identify the elements where the 
interface occurs and correct the appropriate matrix elements in a way described below to 
take into account the effects of the interfaces. Remarkably, these corrections are solutions of 
linear algebraic systems of equations. Thus, the matrix remains sparse and computational 
time is kept to a minimum. 



The original formulation [15J] of the IIM does not consider eigenvalue problems such as the 
problems of interest here. Hence, in order to deal with waveguide problems for Bragg fibers, 
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the method must be extended to handle any eigenvalue problem described by a second 
order differential operator. Moreover, the method must be extended to handle coupled 
equations like the ones describing the two polarization components of the electromagnetic 
field. Furthermore, since the method is based on finite differences one can use higher order 
schemes to increase the accuracy of the calculations. The extension to higher order accuracy 
is also presented in this article. However, extending to higher order posses a major limitation. 
In some problems, the geometry of the interfaces are such that in order to have enough 
points between them (using a uniform grid) requires to increase the total number of points 
and as such computational time. To overcome this, we introduce a coordinate stretching 
transformation which allows the method to handle these more demanding geometries. 

The article is organized as follows: We begin with the description of the method in 
second order. In so doing, we extend the original formulation of the IIM to matrix eigenvalue 
problems. In addition, it is shown that all discontinuities/singularities are removed from the 
equation and passed on the FD scheme as corrections to the standard FD coefficients based 
on matching conditions across an interface. These corrections are calculated using linear 
systems of algebraic equations. Then fibers with deformations are considered to further 
illustrate the versatility of the method. Finally, we extend to fourth order and conclude 
with more demanding geometries in which the original IIM would fail unless a coordinate 
stretching is applied. 



1. Formulation 
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The vector Helmholtz equations in cylindrical coordinates for the magnetic field are [25 
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where k is the wavenumber, n = n(r) is the (arbitrary) index of refraction, and is the 
propagation constant. We focus on the first two equations, since the components of the 
electrical field, E r , Eg and E z , as well as the transverse magnetic field, E z can be recovered 
from E r and Eg using Maxwell's equations. All fields of the guiding modes are assumed to 
go to zero as r — >■ oo. In addition, in order for system (TjQ) to be well defined at the origin, 
the following boundary conditions must hold at r = 
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Separation of variables in Eqs. ([T]) suggests that the fields can be expressed in the form 



EJr, 



E rr 



r cos [rat 



and Eg (r, i 



Eg m (r) sin(m#) with m an integer. Hence, 



Eqs. ([TJ) become 
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In the absence of angular dependence, i.e., m = 0, these equations uncouple and describe 
the TE and TM modes of the fiber, respectively. Hereafter we drop the double subscript 
notation and we set H rm = H r (r) and Hg m = Hg(r). The boundary conditions as r — > oo 
remain the same and at r = Eqs. ([2]) become 

(1 + m 2 )H r + 2mH e = (4a) 
2mH r + (1 + m 2 )H e = (4b) 

When m^l these simply imply that H r (r = 0) = Hg(r = 0) = 0. When m=l, however, 
Eqs. (j3J) are identical and an additional boundary condition must be imposed. It is straight- 
forward to show via a Taylor series expansion around r = that the solution will satisfy the 
condition 

^-(H r -H e ) = 
dr 

which we will impose as our additional boundary condition. The way to implement these 
boundary conditions into the FD scheme is shown in the appendix. 



2. The second order method 

Consider the system of coupled equations that describe the electric, H r , and magnetic, 
H e , fields in a circular waveguide, Eqs. (J3J). Expanding all derivatives in Eqs. ([3]) and after 
appropriate simplifications the system is written as 

d 2 H r ldH r 2m TT f t22 m 2 + 1\ _ o2 TT , r . 
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and the prime (') denotes differentiation with respect to r. Consider a finite difference 
approximation for Eqs. (J5J) of the form (central differences) 

7i# ri i_i + 72#r,i + lzH r ,i+\ + &H d)i = (3 2 H rti (6a) 
SiHe^ + S 2 H e>i + S 3 H 6ji+1 + TH r , = [5 2 H e , (6b) 



where (see appendix) 

_ 1 1_ _ fc2 2 _ J_ _ ^ J_ _J_ _ A _2m 
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with r G [a, b] defined on a uniform grid as 

fb-a\ 

Ti — a + ih — a + il I, z = 0, 1, 2, AT. 

Using this one can find the correct row where the correction must be applied. If the interface 
occurs at r = r* then setting 



j = int 



r — a 



h 



n(r) 



gives the row to be corrected. The function int{} denotes integer part. 

In Eq. (J5]) the index of refraction is a discontinuous function and changes, say, at r = r*, 
so that 

Tlx, r < r* 
n 2 , r > r* 

The above scheme cannot be used without any modifications as it does not take into consid- 
eration the singularities appearing in the equations due to the form of the index of refraction. 
Thus we reformulate the problem, including the differential equations, in vector form. In 
addition, the terms including derivatives of discontinuous functions are neglected (we as- 
sume the index is piecewise constant) and their contribution is incorporated into the finite 
difference scheme through appropriate jump conditions on the interfaces. In so doing, Eqs. 
OS]) read in vector form 

H rr + -H r + BH = (3 2 H (7) 

r 

where H = (H r , Hg) T and 

k 2 n 2 — (m 2 + l)/r 2 — 2m/r 2 

—2m/r 2 k 2 n 2 — (m 2 + l)/r 2 

In the vector formulation subscripts denote differentiation. 

The goal is to determine the coefficients of the finite difference approximation to take 
into account this jump of the refractive index at r = r*. To do this, divide the region [a, b] 
into two, the (— ) region for r < r* an the (+) region for r > r*, as in Fig. [U The analysis is 
similar for the two regions, but needs to be repeated for both. Start with the (— ) region: we 
need to replace Hj_i, H i; Hj +1 in Eqs. (Q so that the local truncation error is first order. 
Expand H around the points before and after the jump, namely 

H(r i _ 1 ) = H i _ 1 = H- + (r H -r*)H r -^(r H -rfH; (8a) 
H(r,) = H, = H - + (r i -r*)H r - + ^(r,--r*) 2 H- (8b) 
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Figure 1: The (— ) and (+) regions, the problematic point r — r* and the irregular grid points rj and rj + \. 



The index j denotes points closest to the jump, as in Fig. [TJ Notice that we only include 
second order terms in the expansions. We need to replace the (+) functions in Eq. flHcj) 
since we are in the (— ) region; this is done through the continuity conditions. 

To derive these continuity or matching conditions one needs to refer to the physical 
properties of the problem. Here all fields are continuous functions across all interfaces, i.e. 
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The last equation and Eq. ffTcj) also yield (recall that all fields are only functions of r) 
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and since all fields are continuous across r = r* the only nonzero remaining terms are 
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Finally, the matching conditions for the second derivatives are a consequence of the conti- 
nuity of the fields and Eqs. ([5]) since 



which results in 
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In summary, the continuity conditions in vector form are 

H = H + 

= CHZ + DH 
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To put everything together return to the FD approximation in matrix form 

riHj_i + r 2 Hj + r 3 Hj + i = /? 2 H; 

where the scalar coefficients 7's are replaced by 2 x 2 matrices. Replacing H 3 -_i, Hj and 
Hj+i from Eqs. using Eqs. (EJ) and 

/3 2 H, = H" + -^H r ~ + B H 

we obtain an equality with H~ , H~ and H on both sides. Matching the relative coefficients 
results in the following linear system for the coefficients at rj < r* 
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For the latter system (the (+) side, r J+ i > r*) Eqs. (Q were inverted to substitute for the 
(— ) side. Each of the above systems represents a 12 x 12 system of algebraic equations that 
determines the coefficients of the matrices. If multiple interfaces are present, one merely 
applies these difference formulas multiple times. Note that the result is a system of finite 
difference equations each involving three neighboring points making the resulting equations 
tridiagonal. Because of the tridiagonal structure of the matrix, sparse matrix algebra can be 
used to determine the eigenvalues and eigenmodes. Thus, a large number of points can be 
used for modest computational cost, which allows the accuracy of the results to be increased 
and the modes of complicated structures to be determined, e.g., Bragg fibers with many 
thin layers j^j. Also note that the corrections depend only on the values of the refractive 
index before and after the discontinuity. This means that the jump conditions do not have 
to be modified if the index varies radially between the discontinuities. 

To test our method we use a Bragg fiber with an air core of radius 1.0 /im and a cladding 
that consists of alternating layers with refractive indices n\ = 3.0 and n 2 = 1.5 0, [l2[. The 
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an imaginary 



distance between layers is 0.130 /im and 0.265 /im. As done previously 
cladding with a refractive index close to zero is added outside the core to prevent reflections. 
For multilayer Bragg fibers the effective index, (3/k, is usually measured instead of just (3. 
Within the spectral range of 1.4 /im < A < 1.6 /im, the Bragg fiber supports a single TE 
mode, whose propagation constant effective index is plotted in Fig. |2j 

This method can also be applied to more complicated and computationally demanding 
fibers. For example, let's consider the Omniguide fiber described in Ref. [5J. This is 
a large air-core fiber with core radius 13.02/iin surrounded by 17 layers, starting with a 
high-index layer, with indices n\ = 4.6 and n 2 = 1.6, and thicknesses l\ = 0.09548/im 
and l 2 = 0.33852/im, respectively. Under these parameters the wavelength for the lowest 
dissipation losses is A = 1.55/im. The first two modes of the fiber, namely the TE i and 
TEu are plotted in Fig. [3J The effective indices are 0.99736632 and 1.00097643, respectively. 
The way to implement the additional boundary condition for the TEu mode (m = 1) is 
described in the appendix. 
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Figure 3: The fundamental TE mode and the TE\\ mode of the air-core Omniguide fiber described in the 
text and in Ref. Q. The bottom figure is a blow-up of the fields near the end of the fiber's core. 



10 



Notably, plastic optical fibers (POFs) which attracted recent attention because of their 
use in subscriber line systems and home networks also have a large core diameter and a high 
core-cladding refractive-index difference compared with conventional silica glass multimode 
optical fibers and can support tens of thousands to hundreds of thousands of propagation 
modes. A recent finite element method was used 0, 0| to analyze their properties. We 
discuss these (and more demanding geometries) at a later section. 



3. Fibers with deformations 

Single defects surrounded by Bragg reflectors as the basis for annular resonators were 
proposed and analyzed in Ref. j26|. The basic geometry was a circumferentially-guiding 
defect is located within a medium which consists of annular Bragg layers. As a result of 
the circular geometry, the layer widths, unlike in rectangular geometry, are not constant, 
and the task is to determine the widths that lead to maximum confinement in the defect. 



In addition, it has been suggested [27J that fibers with such defects can be used to model 
pairs of identical touching hollow Bragg fibers. The dielectric profile along the interfiber 
center line resembles a one-dimensional Bragg grating with a central defect formed by the 
two external layers of the fiber mirrors. 

Figure H] depicts the magnetic field inside a defect. The high index layers and the defect 
have an effective refractive index m — 2 while the low index layers have an effective refractive 
index n 2 = 1. The internal and external Bragg reflectors have 10 periods, and the wavelength 
is 1.45/im. The defect is (A/2)/xm wide. The effective index is found to be 0.92257830. 




rGim) 



Figure 4: The magnetic field distribution of an annular defect mode resonator. 



Resonant features correspond to the points of accidental degeneracy of TEoi with higher- 
order modes. This is true for this case and the effective indices for the TEqi and TE20 modes 
where found to be 0.87238117 and 0.49188825, respectively For more details on the physical 
properties of these defects we refer the reader to [27| . 



11 



4. Higher order method 

Higher order accuracy methods can also be derived in a similar manner. The main focus 
is again Eqs. As before, assume the refractive index to be piece- wise constant, and 

write these equations in vector form as 

H rr + -H r + BiH = (3 2 H (10) 
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and the subscripts denote differentiation. Recall, the continuity conditions are 
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The fourth order finite difference approximation in matrix form is 
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where (see appendix) 
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at all regular points. We need to define these matrices for the irregular points using the 
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IIM. We define our grid as before, namely r« = a + iKf, with N the total number of points, 
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and assume that the point r = r* is between two points, say rj and rj + \. As in the case of 
second order accuracy, we define a point to be regular if all points of the finite difference 
equation, Eq. (Tl2|) . are on the same region, either the (— ) or the (+) regions. All other 
points are irregular. Thus, in Fig. [T]the irregular points are rj_i,rj,rj +1 and rj +2 . We only 
need to define the T's in Eq. (I12p at these points. 

To do so we need to expand in Taylor series all function around the problematic point 
r = r* up to and including terms of fourth order. Hence, for example at r = rj < r* 
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While the continuity conditions are known up to the second derivative, see Eqs. (|lip . 
additional are needed for the higher derivatives. These are obtained from differentiating Eq. 
(JTj) and using Eqs. (TTTI) . Hence for the third derivative 

H^ rr H — + B r H^ = /3 2 H^!~ = /3 2 (C*iH r + DiH) 

which after further use of Eqs. (j7j) and (ITT1) becomes 

H^!~ rr = C^H",. + Z^3H rr . + E 3 H r + -F3H 
and similarly for the fourth derivative 

H = HZ. + C 5 H~ + -DsH" + -E^H" + F 5 H 



where 

^ = ( n 2 /n 2 ; ' 

n ( 
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5 ~ V 2m(n 2 /n 2 -l)/r* 3 2fc 2 (n 2 /n 2 - l)(n 2 2 - n\)/r* - {Am 2 + 10)(n 2 /n 2 - l)/r* 3 
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The elements of the other two matrices 
[Fs]n = ^(nl-nt) + 2^'^Mf- 



are 
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[F 5 ] 12 = -12m^^ 



[F 5 ] 22 = k\nl-n\Y + 2 -^>^ 

Since higher derivative are present we must complete the system of equations by expand- 
ing the right hand side as follows 

/3 2 h, = ^(^ + {r j -T*)n-+ l -{r j -T*fn-;)j 
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where we have differentiated Eq. (fTUj) twice and used the appropriate continuity conditions. 
Also, 
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Now, we can form a system of algebraic equations to determine all the unknown coefficients. 
This is done as before, by equating the coefficients of H~ rrr , H~ rr , H^ r , and H at r < r* 
and their (+) counterparts at r > r*. Thus at r = rj < r* 

r\ + r 2 + r 3 + + i)r 4 + s^j + 2)r 5 = b 1 + - r*)B 3 + ^( rj - r*) 2 B 6 
( r ._ 2 _ r * )ri + ( r ._ : _ r *)r 2 + fa - r*)r 3 + s 2 ( 3 + i)r 4 + s 2 (j + 2)r 5 = 

^/ 2 + (r J -r*)5 2 + i(r,-r*) 2 5 5 
±fa_ 2 - r*) 2 r x + ±fa_! - r*fT 2 + ±fa - r*) 2 r 3 + S 3 (j + l)r 4 + S 3 (j + 2)r 5 = 

±fa_ 2 - r*) 3 r\ + ifa-i - r*) 3 r 2 + i( r ,- - r*) 3 r 3 + S 4 (j + i)r 4 + 5 4 (j + 2)r 5 = 

(r,-r*)/ 2 + ^(r,-r*) 2 / 2 
l(r,_ 2 - r*) 4 r x + l(r,_! - r*) 4 r 2 + Ifa - r*) 4 r 3 + 5 5 (j + l)r 4 + 5 5 (j + 2)r 5 = 

ifa-r*) 2 / 2 

where ^(j) = J 2 + fa— r*)£>i + ±fa - r*) 2 iq + |fa - r *) 3 F 3 + ^fa - r*) 4 F 4 , S 2 (j) = fa- 
r*)C 1+ §fa - r*) 2 £ 1+ ±fa - r') 3 ifc+ifa - r*) 4 £ 5 , 5 3 (j) = Jfa - r*) 2 / 2 +|fa - r*) 3 D 3 + 
^fa - r*) 4 ^, 5 4 (j) = ±fa - r*) 3 C 3 + ^fa - r*) 4 C 5 and 5 5 (j) = £fa - r*) 4 / 2 . At r = 
r j+1 > r* 

S e (j - 1)^ + 5 6 (j)r 2 + T 3 + r 4 + r 5 = B 1 + fa +1 - r*)B 3 + ^fa +1 - r*) 2 5 6 
SrU - m + 5 7 (j)r 2 + (r j+1 - r*)T 3 + (r, +2 - r*)r 4 + fa +3 - r*)r 5 = 
^/ 2 + (r, +1 -r*) J B 2 + i(r, +1 -r*) 2 J B 5 

5 8 (j - l)r\ + S 8 (j)T 2 + l(r j+1 - r*) 2 r 3 + l -{r j+2 - r*) 2 r 4 + ifa +3 - r*) 2 r 5 = 

S 9 (j - l)r\ + S 9 (j)T 2 + i(r, +1 - r*) 3 r 3 + ifa +2 - r*) 3 r 4 + ifa +3 - r*) 3 r 5 = 

fa +1 -r*)/ 2 + ^fa +1 -r*) 2 / 2 
Sioti ~ m + S w (j)r 2 + ^fa + i - r*) 4 r 3 + l(r i+2 - r*) 4 r 4 + l(r i+3 - r*) 4 r 5 = 

^fa +1 -r*) 2 / 2 
15 



where S e (j) = I 2 + (r 3 -r*)D 2 + 



F2 + \{r ] -r*fF, + l 



24 V r i 



,*\4 



Fe, S 7 (j) 



r*)C 2 +\{r 3 - r *) 2 E 2 +\{r 3 - r*) 3 E 4 +±( rj - r*) 4 £ 6 , S 8 (j) = \{r 3 - r*) 2 I 2 +\{r 3 - r*fD A + 



6 ^* 3 

L( rj -r*)±D 6 , S 9 (j) = l(r, 



6' 



r 



j C 4 + ^^j- 



C 6 and Si (j) 



24 l r J 



,*\4 



J 2 . The 



_i and r = r 7 _ 2 are obtained in a similar manner and will not 



systems at the points r 
be explicitly given here. 

To demonstrate the convergence of the method we choose a fiber with multiple layers 
that can support both TE and TM modes. Consider a fiber consisting of an air core of 
radius 5/mi and a series of alternating layers of radii 1/im and refractive indices n\ = 2 and 
n 2 = 1, respectively. The results for the propagating constant for the TE and TM modes are 
shown in Table [TJ In Fig. the corresponding TE and TM modes are shown. The m = 1 



Mode 



Propagation constant 



N=500 



N=1000 



N=2000 



N=4000 N=10000 N=20000 



TM 
4th order 

TE 
4th order 



0.52520655 
0.52732576 
0.78582132 
0.78334621 



0.52731320 
0.52707445 
0.78580424 
0.78575156 



0.52714039 
0.52705720 
0.78588272 
0.78588122 



0.52707988 
0.52705650 
0.78590279 
0.78590670 



0.52706043 
0.52705645 
0.78590848 
0.78590941 



0.52705446 
0.52705645 
0.78590929 
0.78590954 



Table 1: The propagation constants for the multi-layered fiber of the text. 



case is again described in the appendix. 




5 10 15 20 25 30 35 40 

r(nm) 



Figure 5: The TE and TM modes respectively corresponding to the propagation constants of Table [TJ 



One can go to even higher higher accuracy. In fact in Ref. [28J the authors have obtained 
sixteenth order accuracy. This, however, introduces a fundamental limitation, namely the 
number of irregular points. In order for the method to work, one must have a minimum 
number of points between interfaces. When the structure of the fiber become more compli- 
cated it is obvious that this is a serious limitation. We demonstrate how to overcome this 
in the following section. 
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5. Coordinate stretching transformation 



P 



One of the problems one might encounter in this formalism is the total number of points 
in a grid that must be used to have high accuracy results. The number of grid points is 
directly proportional to the size of the matrix to be diagonalized, thus more points means 
more computational time. When the width of the fiber's core is considerably (orders of 
magnitude) larger than the width of the layers, such as in the case of the fibers in Refs. 
[I- S, 0, 0, 0] , one needs to use a sufficiently small numerical step size to insure that there 
exist at least two points per layer. If the number of points per interface is less than two, the 
assumption used in the derivation of the method is not true, and the method fails. This is 
because the method assumes that there is only one interface for each set of three grid points. 
Because the reflecting layers are small, this means that a small step size must be used. If 
one uses a uniform mesh, this means that in the core region one is using a step size that is 
much smaller than necessary, i.e., computational effort is being made unnecessarily. 

To avoid this problem, we introduce a coordinate stretching by using a new independent 
variable p such that 

r, r < R* 
R* + a(r-R*), r > R* 

Here R* is the location of an additional artificial layer placed arbitrarily inside the core of 
the fiber and a is a stretching parameter such that a > 1. This transformation increases 
the effective width between the layers relative to the core (by a factor of cr). We will use 
uniformly spaced points in terms of p, which is equivalent to using a smaller step size in 
the core and a larger step size in the layers when measured in terms of r. Under this 
transformation, Eqs. ([7]) become 

H pp + -H„ + BH = /3 2 H, p<R* (13a) 
P 

2 

a 2 H pp + ° U P + BU = (3 2 U, p > R* (13b) 

p + (a - ljR* 

At the point p = R* we impose the additional artificial jump conditions 

H + = H , = H r 

so that, in terms of the new coordinate, both fields are continuous and their derivatives 
satisfy 

ctH p = H p 

It follows then, from Eq. ([7]) and the above jump conditions that the second derivatives 
satisfy 

where the (— ) and (+) regions are on the left and right of the discontinuity, respectively. 
At the real layers, including the core, Eqs. ([9]) apply, with appropriate changes, namely 

P-R* , y->* d d 

r = + R, — = a— 

a dr dp 
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The finite difference approximation in matrix form at all points is again 

TiHj.i + r 2 Hj + r 3 H i+ i = /3 2 Hj 

where the T's are 2x2 matrices. As always, we need to expand Hj_i, Hj and Hj+i around 
the points of the grid that include the discontinuity. Assuming that r\, < R* < r J+1 , the 
systems that determine the finite difference coefficients at these regions then satisfy the 
following systems 
At Pj < R* 

T 1 + T 2 + T 3 = B 



- R*)^ + ( Pj - R*)T 2 + r 3 (p i+1 - R*)/a = —I 2 



^{Pi-x - Wfv x + ~(p, - ir) 2 r 2 + ^(p j+ i - R*) 2 r 3 = h 



At p j+1 > R* 



T 1 + T 2 + T 3 = B 
a( Pj - i?*)r x + (p j+1 - R*)T 2 + T 3 (p J+2 - R 



a 
R* 

*\2 r _2 ; 



-iPi - wyr, + -( P , +1 - R*yr 2 + -( Pj+2 - R*yr 3 = *% 

The finite difference coefficients for the rest points are determined as usual to be, at pj < p* 



Ti + r 2 + r 3 



h + {Pm- P*)D+ l -{p 3+1 - P *) 2 F 



B~ 



{p J _ 1 -p*)T l + {p J -p*)T 2 + T 3 



1 



( Pj+1 -p*)c+~(p j+1 - P *yE 



-( Pi _x - p*fY x + ~( Pj - p*yT 2 + ~(p j+1 - p*YT 3 = a 2 L 



and pj > p* 



I 2 -( Pj -p*)C- l D + -{p 3 -p*yF 2 



+ r 2 + r 3 = b + 



( Pj -p*)c- 1 + -( Pj -p*) 2 j e; 2 



+ (pj+i - p*) v 2 + r 3 (p i+2 - p*) 



1 



i 



i 



{ Pj - p*yr x + -(p i+1 - p*)^r 2 + -( Pj+2 - P *yr 3 = 0*1 



a 



P* 



2 ' ' 2 ' ' 2 

where the matrices E and F are defined in section [2] and we need to introduce the matrices 

E2 = ' (1 - n\/nl)/p* ) ' F2 



P(n2-n 2 ) 
m{\ - n\/nl) / p* k 2 {n\ - n{) - (1 - n\/n 2 2 ) / p* 2 
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Using the coordinate stretching more demanding structures can be analyzed using the 



same number of points. Consider, for example, the fiber described in Refs. |29l. l30l] . This 
fiber consists of a hollow-core 316 fim in radius surrounded by 70 alternating layers of As2Se3 
0.27 /zm thick, and polyether imide 0.47 /im thick. The fundamental photonic bandgap is 
centered at A = 2.28 /im and the refractive indices of the layers are n\ = 2.8 and n 2 = 
1.55 + 1. Ox 10" 4 i, respectively. The effective indices are found to be 0.99999032 + 1.4xl0~ 12 i 
and 0.99999025 + 7.1xl0~ n i for the TE and TM modes, respectively. These results were 
obtained for 20000 points. Now consider the coordinate stretching transformation: 



P 



r, r < 300/zm 
300/zm + 5(r — 300/xm), r > 300/zm 



For 5000 points the IIM in its original formulation fails since the three points per interface 
fails. However, we can get results of good accuracy using the above transformation. In- 
deed, for 5000 points the results are for the effective index 0.99999030 + 1.39 x 10 _12 i and 
0.99999021 + 6.9 x 10 11 i for the TE and TM modes, respectively. In Fig. [6]one can see and 
compare the TE modes obtained with and without the coordinate stretching. It is apparent 
that part from the added interface the fields are identical. This transformation proves very 
useful for demanding structures such as this one. 



6. Conclusions 

We presented a numerical method based on the immersed interface method that can be 
used to obtain the propagation modes of circularly symmetric Bragg fibers with arbitrary 
index profiles. In its original formulation the method is second order accurate and was 
applied to boundary value problems with discontinuous and/or singular coefficients. We 
extended this method to matrix eigenvalue problems and to higher accuracy. Cumbersome 
integrations or finding roots of nontrivial functions, such as Bessel functions, are avoided and 
computational time is minimized without sacrificing accuracy. All modes can be determined 
and excellent results are achieved even for fibers with complicated structure. Even when the 
geometry of the fiber is rather demanding, as in Omniguide fibers, a coordinate stretching 
can be applied to keep computation time to a minimum. 
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Appendix A. Finite difference approximations for first and second derivatives 

For the second order scheme 

dH H i+ i — ifj-i 

~dr ~ 2h 
d 2 H H i+ i — 2Hi + Hi-\ 

dr 2 h 2 
19 
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Figure 6: The TE mode of an Omniguidc fiber using the original IIM and the coordinate stretching trans- 
formation. In the first figure only the first and last layers are plotted. 
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and for the fourth order case 

dH —H i+ 2 + 8H i+ i — 8i/j_i + Hi-% 
~~d~r~ ~ 12h 

d 2 H -H i+2 + lRgj+i - 30ff, + 16^-i - H^ 2 
dr 2 12h 2 

Appendix B. Boundary conditions for higher-order differences 

As usual, special care must be given to the end points. High-order methods require an 
extrapolation scheme to determine the differencing coefficients at r 2 and r^-i- Given the 
boundary conditions H(ri) = Hi = and H(r N ) = H N = we use the extrapolation 

H 2 = -H 6 + 4H 5 - 6H 4 + 4H 3 

Hjv_i = 4Hiv_ 2 — 6Hjv_ 3 + 4Ha?_ 4 — Hns 

to accommodate for these extra points. 



Appendix C. Special case boundary conditions (m = 1) 

All of the cases considered had the fields tend to zero at the boundaries, namely at r = 
and r — > oo. However, this is not the case for the calculation of the HEu mode (our m = 1 
case). Indeed, recall that the boundary conditions for this case are 

H r + H e = and ^- (H r - H e ) = 

dr 

at r = 0. To accommodate this we need to alter the finite difference scheme to include the 
new boundary conditions. 

Appendix C.l. The second order correction 

At the first point where i — 1 Eqs. ([6]) become 

7l#r,0 + 72#r,l + 73#r,2 + A# e>1 = ^E r ^ (C.la) 

5 x He tG + 5 2 H 0A + 5 3 H 6 , 2 + VH r>1 = (5 2 H e ^ (C.lb) 

and the boundary conditions become under a second order approximation of the first deriva- 
tive 

H r ,i + He,i = and H r>2 — H r>0 = H Q2 — Hg :0 

Note that all the coefficients are known (we are not at an interface) and we only need to 
eliminate the two H r fl and ifg.o terms. 

Adding Eqs. ( 10. ip and taking the limit as r — > gives (recall 71 = Si, 72 = S 2 , 73 = S 3 
and T = A as clearly seen from their definition below Eqs. flS}) 

H r2 + H r0 = H d2 + Hqq 

which suggests that 

H = H 2 
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Appendix C.2. The fourth order correction 

At the first point where % = 1 we write the system as 



7l#r,-l + 72#r,0 + 73#r,l + 74#r,2 + 75#r,3 + A# M = /3 2 # r ,l (C.2a) 
+ 5 2 #0,O + <f 3 #<M + <M*«,2 + 4^,3 + r# ril = /3 2 ^,! (C.2b) 

where the definitions for the coefficients follow from Eq. ffl~2]) . The steps here follow the 
procedure for the second order results. Thus, from the definition of the first derivative at 
fourth order for i = 1 we get 

H r -\ — 8H rt0 + 8H r2 — H rS = H e _i — 8H 9 fi + 8H 6:2 — Hq i3 

As above adding Eqs. (10. 2[) and taking the limit r — > gives 

^(flr.-i + #0,-0 - 2 -{H Tfl + F 0)O ) + ^(# r , 2 + H d}2 ) - ^(H r fi + H 9 fi) = (C.3) 

However, we also need to extrapolate for the point i = — 1, namely 

H_x = 4H - 6Hi + 4H 2 - H 3 (C.4) 

Note that while the point % = was before ignored this is not the case here. Using these 
equations we can now eliminate the point i = as follows 

3 1 
H = — -Hi + 3H 2 — -H 3 
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